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Abstract 



We extend previous results (arXiv: 0804.2630 [hep-ph]) on factoriza- 
tion in high-energy nucleus-nucleus collisions by computing the inclusive 
multigluon spectrum to next-to-leading order. The factorization formula 
is strictly valid for multigluon emission in a slice of rapidity of width 
C^^ ' AY < aj^. Our results shows that often neglected disconnected graphs 

^^ , dominate the inclusive multigluon spectrum, and are crucial in order to 

achieve factorization for this quantity. These results provide a dynamical 
framework for the Glasma flux tube picture of the striking "ridge" -like 
correlation seen in heavy ion collisions. 



^ '. 1 Introduction 



In a recent work, henceforth referred to as Paper I [1], we investigated the formal 
basis for the application of the Color Glass Condensate (CGC) framework [2-8] 
to the collision of two high energy nuclei. In Paper I, we focused on the for- 
malism to compute the single gluon inclusive spectrum in the Leading Log x 
approximation. The main result of Paper I is a proof that terms containing lead- 
ing logarithms of l/a;i.2 that arise in all order loop corrections to this spectrum 
can be factorized in the distributions of color sources T4/^[pi.2] in each of the two 
nuclei, evolved with the JIMWLK equation [9-16] from the beam rapidity to 
the rapidity of the measured gluon. One obtains for the single inclusive gluon 



distribution the result 

dN 



^) = J[Dpi][Dp2]Wy,^^^-y[PiWy,..^+y[P2] ^ 



(1) 



The W functionals are universal properties of the nuclear wavefunctions at high 
energies and (in analogy to the parton distribution functions of coUinear factor- 
ization) can be extracted from deep inelastic scattering or proton-nucleus scat- 
tering experiments off nuclei. The inclusive single gluon spectrum (dN/cfip) 
that appears under the integral in the right hand side is the Leading Order spec- 
trum corresponding to one configuration of the sources pi^2 ~ it is obtained by 
solving the classical Yang-Mills equations for this fixed distribution of sources. 
This factorization theorem allows for considerable predictive power by relating 
measurements in a variety of scattering processes. It should be particularly 
useful at the LHC, where the rapidity reach in proton- nucleus and nucleus- 
nucleus collisions will be considerable and the effects of energy evolution of the 
distribution of color sources clearly visible. 

The derivation of the factorizcd expression in eq. (1) relied on two essential 
steps: 

1. The 1-loop corrections to the gluon spectrum can, in the leading logarithm 
approximation, be expressed as the action of a certain linear operator on 
the leading order spectrum^ , 

2. This operator acting on the initial color fields on the light-cone is, again 
in the leading log approximation, the JIMWLK Hamiltonian. 

In the present paper, we will show that a straightforward generalization of the 
first of these two steps is sufficient to extend our factorization result to inclusive 
multigluon spectra when all the measured gluons are located in a rapidity region 
of maximal width Al" < a~^ . 

The paper is organized as follows. In section 2, we define a generating func- 
tional for multiparticle production in nucleus-nucleus collisions. This extends 
to the QCD case our previous results [17,18] for a similar object introduced for 
a (JT' theory. We discuss key features of this generating functional and develop 
a diagrammatic interpretation of this object. We show how (at leading order) 
its first derivative can be expressed in terms of classical solutions of the Yang- 
Mills equations that obey both advanced and retarded boundary conditions. In 
section 3, we consider in detail the inclusive two-gluon spectrum. We obtain an 
expression of this spectrum at next-to-leading order (NLO) using the previously 
defined generating functional. We end the section by showing that the leading 
logs of l/a;i^2 in this quantity can be factorizcd in the distributions of incoming 
color sources, provided the rapidity separation between the two gluons is small 
enough. We show that our formalism gives rise to the Glasma flux tube pic- 
ture [19], which has been suggested as a mechanism to describe the ridge-like 
structure observed in heavy ion collisions at RHIC [20-24]. In section 4, we 



^See the discussion after cqs. (40-41) and at the end of section 3.5 in [1]. 



generalize this factorization result to the case of the inclusive n-gluon spectrum. 
Knowing all the moments defines the complete probability distribution. We 
demonstrate how the leading logarithmic corrections to the multiplicity distri- 
bution can be factorized into the JIMWLK evolution of the sources. We end 
with a brief summary. The three appendices are devoted to the more technical 
aspects of our discussion. 

2 Generating functional 

In Paper I, we developed the tools for studying at LO and NLO the single in- 
clusive gluon spectrum in AA collisions in the CGC framework. Our goal is to 
generalize these techniques to obtain similar results for the n-gluon spectrum. 
Towards that purpose, we will define in this section a generating functional for 
n-gluon production, discuss its properties and develop a diagrammatic interpre- 
tation. We then discuss the LO computation of the first derivative of this object 
in terms of solutions of classical Yang-Mills equations with both retarded and 
advanced boundary conditions. 

2.1 Definition and properties 

We define the generating functional as 

J^^ 'iff 2 

where we use the following compact notation for phase-space integrals^, 

f ^ f d^P .3. 

J -J (2^)32i?p • ^""^ 

p 

In this definition, z{p) is an arbitrary function over the 1-gluon phase-space. 
The matrix clement squared that appears in the right hand side is implicitly 
summed over the polarizations and colors of the produced gluons. Note that in 
this section, we consider the external current J^ coupled to the gauge field to 
be fixed. This is the case in the Color Glass Condensate (CGC) framework [8] 
where the fixed sources represent the large x light cone color charge densities 
in the nuclear wavefunctions. We will address the issue of averaging over the 
external color sources later in this paper. 

The generating functional generalizes the generating function F{z) we intro- 
duced in ref. [17]. This previously defined function is simply obtained as 

T[z{p) ^ z*] = F{z*) , (4) 



■^Whenever the integrand contains po in such integrals, it should be replaced by the positive 
on-shell energy po = \p\- 



when z{p) is a constant z* . Another obvious property of JF[z(p)] is 



T[z{p) ^ 1] = 1 



(5) 



which is a consequence of the fact that the theory is unitary. 

The generating functional encapsulates the entire information content of the 
nuclear collision within the CGC framework. Indeed, if JF[z(p)] were known, 
one could use it to build an event generator for the early Glasma [25,18,8] 
stage of nucleus-nucleus collisions. In particular, one can compute the inclu- 
sive multigluon spectra. For instance, the single inclusive"^ gluon spectrum is 

obtained as 

dN __ ST[z] 

(Pp 



Sz{p) 



(6) 



Likewise, the inclusive 2-gluon spectrum is obtained by differentiating J-'[z] 
twice. 



dm 



S^T[z 



d^pd^q Sz{p)Sz{q) ^^^ 

where the integral over p and q on the left hand side of this expression is the 
average value of N{N — 1). Physically, this quantity, in an event, corresponds 
to a histogram of all pairs of distinct gluons with momenta (p, q) . We will 
discuss the average over all such events later. Eqs. (6) and (7) are the two 
simplest examples of the use of this generating functional, but in principle one 
can derive from it any observable that is related to the distribution of gluons 
produced in the collision. Eq. (7) can be generalized to 



d'^Nr, 



d"T[z] 



d^p^ ■ ■ ■ d^p^ 5z{p-^) ■ ■ ■ dz{pj 



(8) 



for the inclusive n-gluon spectrum. Note that the l.h.s, integrated over the n- 
particlc phase space, is normalized to the average value of N(N — 1) • • • {N — 
n+l). 

From eq. (8), it is possible to represent the generating functional J-'[z] as 



oo _. -. n 



(i"iV„ 



d^p-^ ■ ■ ■ d^p„ 



(9) 



This formula will later be the basis of our strategy to obtain an expression for 
J-"[z] at Leading Log. We will first obtain Leading Log expressions for the n- 
gluon spectra*, and will show that the infinite sum in eq. (9) leads to a very 
simple expression. 



•^Note that setting z{p) to zero instead, after taking the functional derivative, one obtains 
the differential probability for producing exactly one gluon in the collision, 

dPi _ SJ^[z] 

d^p Sz{p) ^=0 ' 



■'With the important limitation that the n gluons all sit in a rapidity slice of width AY < 
-1 ^ 



Once we know T\z\ (with a given accuracy), one can use the fact that its 
Taylor coefhcients at zijp) = are the differential probabilities for producing a 
fixed number of particles^ , 






d^Pu 



d?v^ ■ ■ ■ d^p^ ' 



(10) 



From this second representation of J-'[z], one can extract from J-'[z] detailed 
information about the distribution of produced gluons. 



2.2 Diagrammatic interpretation of J-'[z] 

In order to see what are the diagrams that contribute to J-'[z], let us first define 



V 



Vr. 



(11) 



with 



Vp^j2<(py>^(py Jd'^d^y^^'''^'"'^ °-°y j^, 



x)5r{y) 



(12) 



This operator has already been introduced in [17,18] to write P„ in terms of 
vacuum diagrams. The only difference here is that we extend its definition to 
the case of vector particles and QCD. The sum over the gluon polarizations 
A spans the two physical polarization states. By mimicking the manipulations 
performed for scalar fields, one can prove that 



T[z{p)] ^ exp 



z(p) Vp 



mJt] e-'VlJ':] 



j'^=j'^=.ji^ 



(13) 



where iy[J^] is the sum of the connected vacuum diagrams evaluated with the 
external current J^. It is easy to check that all the formulas we previously 
obtained in [17,18] for P„ or for the generating function F{z) are all particular 
cases of this formula. 

From the interpretation of the operator T) as an operator that makes cuts 
through vacuum diagrams, we see that the functional JF[z(p)] is the sum of all 
the cut vacuum diagrams (connected or not) in which every cut propagator with 
momentum p is weighted by z{p). Let us call iWlJ^, J^; z] the sum of all such 
connected diagrams (before the currents J^ and J'^ are set equal to the physical 
value J^) : 



exp 



z{p) Vp 



^y[J^] ^^iV'[Jt] 



(14) 



^Note that there is no 1/n! in this formula. A quick way to convinee oneself that this is 
correct is to set z{p) = 1; the integrals over the momenta p^ give the total probabilities P„, 
which add up to unity. 



It is useful to compute the first derivative of jF[z(p)] with respect to z{p), 



5J^[z] _ 1 

Sz{p) ^ {2tj:Y2E^ 



jy^W[J-,J-;z] 



(15) 



Performing explicitly the derivatives contained in eq. (12), this can be rewritten 



as 



ST[z] 
Sz{p) 



(2 



^^^ ^ <(p)eX(p)* / d^x d'v e^^<^-y) n.n. 



P A 

5iW SiW 



S^iW 



5J^{x)6J-_{y) 6J^{x)5J'^{y) 



^iW[J'^,Jt;z] 



(16) 



The final exponential in this formula is nothing but T[z\ itself. Therefore, we 
can write 



5\tlT[z] 
Sz{p) 



(2 



-1-—J2e';^{p)el{pr J d^xd^y e^P<-y^ D^Dy 



"P A 

SiW SiW 



S^iW 



SJ^{x) 5J-_{y) 5J^{x)5J-_{y) 



(17) 



This formula tells us that this quantity is made up of only connected diagrams 
since iW is a sum of connected diagrams. We also observe that this formula is 
very similar to the formula for the single inclusive particle spectrum with one 
very important difference: the function z{p) is not set to 1 at the end, and 
therefore appears as a multiplicative factor attached to each cut propagator. 



2.3 5\\iJ^\z\/5z{p) at leading order 

Let us now show that, in the regime of strong external color sources, the ex- 
pression in cq. (17) can be expressed at leading order (LO) in terms of classical 
solutions of the Yang-Mills equations. 

First of all, note that the first derivatives 5W /SJ'^ are of order® g^^, while 
the second derivative S'^W /SJ^SJ'^ is order g^ . Thus the first term, composed 
of the product of two first derivatives, is the leading one. The second term 
begins to contribute only at next-to- leading order (NLO). At LO, we can thus 
write 



S\nT\. 



5z{p) 



{2tt 



;i^^ Yl eAM(p)eA^(p) / d'^xd 



'^y e*^'^"^"^^ 



xD^DyA^+{x)A''_{y), 



(18) 



^Because W is the sum of connected vacuum graphs, in the presence of external sources 



where we denote^ 

tree 5J^{x) 



(19) 



j'^=j'i=ji^ 



The "tree" here means we keep only tree diagrams in the expansion of SW /SJ±. 
that defines ^^. 

AU the arguments developed to compute the generating function F(z) at 
leading order [18] can be extended trivially to the present situation, and one 
obtains the following results : 

• A'^ is a solution of the classical Yang-Mills equations, 

[P^,^^-^] -J% (20) 

• If one decomposes A'^{x) in Fourier modes, 

-*'.« p 

(21) 
with a±p\aix) = e^(p)T"e^*P'^, the boundary conditions obeyed by the 
classical field A'^ (x) can be expressed as simple constraints on the Fourier 
coefficients®, 

/j+\-oo;pAa) == f'r\-oo;pXa) = , 
f[+\+oo;p\a) = z{p) f^+\+oo;pXa) , 
f[-\+oo;pXa) = z{p) f^s\+^-p\a) . (22) 

We see that the dependence of the classical fields A'i^ on the function z{p) 
comes entirely from the boundary conditions^, since the Yang-Mills equations 
themselves do not explicitly contain z{p). In terms of the Fourier coefficients 
fj^ , eq. (18) reads 



5z{p) 



A, a 



"^ Ali^(x) depends on function z(p) as well but we have omitted it from the notation to keep 
notations compact. 

*The derivation of this result is analogous to the scalar case discussed in detail in section 
4.2 of Ref. [17]. 

Note that, when z{p) = 1, the boundary conditions in eqs. (22) become 

/_^+'{-oo;pAa) = /i"'(-oo;pAa) = , 

/i+'(+oo;pAa) = /j+'(-}-oo;pAa) , /j_"^(+oo;pAa) = f^S\+oo; pXa) . 

The two conditions at a;" = +co imply that A+ (x) = A.- {x) everywhere. The two conditions 
at x" = — oo then imply that Wra^o^—^ A±(x) = 0. Therefore, when z(p) = 1, the two 
classical fields A'^ become identical to the retarded classical field with a vanishing initial 
condition in the remote past, and eq. (18) gives the single inclusive gluon spectrum as expected. 



Note that it depends only on the Fourier coefficients of the fields at a;'^ = +00. 
Eqs. (20), (22) and (23) do not provide a practical way to obtain the LO 
generating functional !F[z{p)\ because the solutions depend on boundary condi- 
tions at both ±00. It is not known at present how to solve Yang-Mills equations 
with simultaneous advanced and retarded boundary conditions. Nevertheless, 
the procedure outlined here provides a powerful theoretical tool to compute 
other quantities, that can be obtained as derivatives of the generating func- 
tional. A concrete illustration of this strategy is revealed in the case of the 
2-gluon spectrum in the following section. 

3 Two-gluon inclusive spectrum 

In this section, we will specialize our discussion of the generating functional in 
the previous section to the 2-gluon inclusive spectrum at LO and NLO. We will 
demonstrate that, just as in the case of the single gluon spectrum discussed in 
Paper I, the leading logarithm contributions that arise at NLO can be absorbed 
in the JIMWLK wave functionals of the two nuclei, provided the rapidity sep- 
aration between the two gluons is small enough. As in Paper I, one obtains 
a factorized expression for the leading log 2-gluon inclusive spectrum. In the 
following section, this result will be extended to multigluon spectra. 

3.1 Leading Order 

The inclusive 2-gluon spectrum is obtained by taking the second derivative of 
the generating functional !F[z\^ and by setting the functions z{p) and z{q) to 
unity afterwards (see eq. (7)). Alternately, it is easy to obtain this derivative 
from the derivative of lnJF[z]. We get 



(PN2 _ 5\nT[z]6\nT[z\ 5'^\nT[z] 



(Pp (Pq Sz{p) Sz{q) 5z{p)6z{q) 



(24) 



The first term is simply the product of two single gluon spectra (see eq. (6)), and 
therefore corresponds to the disconnected (independent) production of a gluon of 
momentum p and a gluon of momentum q. In contrast, because In jF[z] contains 
only connected diagrams, the second term corresponds to the two gluons being 
produced in the same graph. Note that these expressions correspond to the 
2-gluon spectrum for a fixed configuration of the external sources pi.2. When 
we average over these sources, some graphs that were disconnected prior to 
averaging become connected. Therefore, even the first term in eq. (24) can lead 
to correlations in the measured 2-gluon spectrum. 

The two terms in this expression do not begin at the same order in g^. In 
our power counting, 

\nT[z] = — Tco + ci g2 + C2 g'^ + • • • 1 . (25) 



This implies that the first term in eq. (24) is of order g~^, while the second term 
is of order g~'^ only. For the 2-gluon spectrum, "leading order" therefore means 
^, and we simply havc"'^'' 



9 



<PN. 



(PpcPq 



dN 
(Pp 



dN 
d^ 



(26) 



No new computations are necessary here because we know how to express the 
single gluon spectrum at LO in terms of classical solutions of the Yang-Mills 
equations with retarded boundary conditions. Note that at this order the —N 
term contributing to N2 is subleading relative to the N'^ contribution because 
it starts only at the order g^^ therefore does not appear on the right hand side 
of eq. (26) which is of order .g^^. 

3.2 Next to Leading Order - I 

We shall now study the inclusive 2-gluon spectrum at NLO~the contribution at 
order g^^ in our power counting. At this order, the tree level contribution to 
second term in eq. (24) must be included. We can write therefore 



d^N2 



d^pd^q 



dN 
d^p 



dN 
d^ 



dN 
d^p 



dN 
d^ 



S'^\nT[z 



5z{p)5z{q) 



(27) 



The first two terms again do not require a new computation because we studied 
in great detail the single gluon spectrum at NLO in Paper I [1]. In particular, 
we recall here the previously derived formula 



dN 
d^p 



d^u[(3-T^ 



+ 9 Zl / / d^^d"^^ [a-kxa ■ T„] [a+kxa ■ T^] 



k E 



dN 
d^p 



£2 



+AiV,,o(p) 



(28) 



In this formula, a±kXa denotes small field fluctuations that propagate over the 
classical field A. The subscripts indicate that these fluctuations begin in the 
remote past as plane waves of momentum ±fc, polarization A and color a. Sim- 
ilarly, /3 is also a small field fluctuation propagating on top of A, but this 
fluctuation has a vanishing initial condition in the past and is driven by a non 
zero source term. S is a surface on which the initial value of the classical fields 
are defined, and d^u is the measure on this surface. The operator Tu is the 



^"One should keep in mind therefore that "LO" corresponds to different powers of g^ for 
the single and double inclusive gluon spectra. 



generator of translations of the initial field at the point it G S. AN-^-^^{p) is a 
term contributing to the full expression. It will not be explicited further because 
it does not contain a leading logarithmic contribution-see the discussion of this 
term in [1]. Because we are interested here in these leading log contributions, 
this term will be dropped in all further equations in this paper. 

At this point, we can rewrite the first two terms of the r.h.s. of eq. (27) as 



dN 
d?p 



dN 



dN 



dN 






dN 


dN 


d^q 

LO ^ 


NLO 


C1+C2 


disc d^p 


d^q 

LO ^ 



(29) 



where the subscript "disc" added to the operator between the square brack- 
ets indicates that when the combination TuTt, in £2 acts on the product 
{dN / d^p){dN I d^ q) , we keep only the terms where the two T's act on the same 
factor^ ^. The subscript here reminds us that these terms are disconnected con- 
tributions that are the product of a function of p and a function of q. 

3.3 Next to Leading Order - II 

The third term of cq. (27), involving the second derivative of the log of the 
generating functional, is new and will be computed here. Fortunately, we need 
this term only at leading order-i.e. 0{g~'^). Therefore, our starting point in 
evaluating this term is cq. (23). Differentiating this equation with respect to 
z{q)^ we obtain 



5'^\nT[z] 



6z{p)5z{q) 



1 



{2tt)32Ep 



X,a 



- Sfi+\+^;pXa) 
5z{q) 



/i \+OQ]p\a) + /|+)(+oo;pAa) 



5f[_ \+oo; pXa) 



5z{q) 



(30) 



Further, differentiating eq. (21) with respect to z{q), one observes that the 
quantities 5fe {+<x>\p\a) / 5z{q) arc the Fourier coefficients of the field 



bUi^) ^ 



_ SAiiix) 



6z{q) 



(31) 



at x" 



hoo. The equation of motion obeyed by this object can be obtained by 
differentiating, with respect to z{q), the equation of motion for A'^. In order 
to do this, it is useful to start from the Yang-Mills equations written in a form 
that separates explicitly the kinetic and interaction terms^^ as 



□.5m' -^r./.^. ^^(^) 



dU{A,) 

dAe^i,{x) 



j:: 



(32) 



''[^2]^,^^AB=[C2A]B + A[C2B]. 

^^Note that the diflferentiation with respect to z(q) does not modify the gauge fixing con- 
dition, provided it is Uncar. Thus, b^^q obeys the same gauge condition as A'^ ■ 
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where U{A) is the Yang-Mills potential in a gauge with a linear gauge condition. 
Differentiating this equation with respect to z{q), we get 

□.,/^a.,a.-- ^^JJ^_;),^^^ ]5^,(.)=0. (33) 

In other words, b'^ „{x) obeys the equation of motion of small fluctuations prop- 
agating on top of the classical field Ae- The boundary conditions necessary 
in order to fully determine 6Jf„(x) are easily obtained by differentiating the 
eqs. (22) with respect to z{q): 

&Vtq(-«^;pAa) = b''Sli-^;p\a) = , 

&L+^(+oo;pAa) = z{p) 6^+^(+oo;pAa) + S{p - q)f[+\+oo;pXa) , 

&^^^(+oo;pAa) = z{p) fe'f ^(+oo;pAa) + S{p - q)f'r\+oo;p\a) , (34) 

where we have introduced the obvious notation 

&<:'^(-";pAa)^^^^^^^ (35) 

for the Fourier coefficients of 5^„. We see that we have non homogeneous bound- 
ary conditions, which will lead to a non zero 5(f despite the fact that this fluc- 
tuation obeys an homogeneous equation of motion. Note also that at this point 
we can safely set z{p) = 1 since we do not need to differentiate with respect to 
z{p) again. This leads to the simplification that when z(p) — 1, the classical 
fields A'^ and A'^ become identical-as can be checked from their boundary con- 
ditions (see footnote 9). In fact, their common value is nothing but the classical 
field that vanishes when x'^ -^ — oo. We will simply denote by A^ the common 
value of these two fields and /(^'(a:;°;pAa) its Fourier coefficients. 

Obviously, eqs. (34) are not simple retarded boundary conditions. Our task 
is now to relate the fluctuations 6^ and their Fourier coefficients to fluctuations 
that satisfy simple retarded boundary conditions. In order to achieve this, let 
us again use the small field fluctuations a'^i^xa- They obey the equation of 
motion (33), and the boundary conditions 



a 



M (^\ _ ^lJ.(i,\rpa^±ik-x 

±kX: 



^{x) = 6^(fc)T"e±*% (36) 



Note that the fields a^i,xa introduced earlier are the analogue of the a'^j^xa 
in the absence of a background field. From this definition, a^k\a has only 
negative energy components at x^ -^ — oo, while a^kXa has only positive energy 
components in this limit. Moreover, the fluctuations a^^^Q provide a complete 
basis for the small field fiuctuations that obey eq. (33). From the boundary 
conditions of a^ at x'^ = — oo, we see that we must have 

bl,ix) = J2j^±'>lkXai^)- (37) 

A, a . 
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The coefficients 7^^" in these hnear decompositions do not depend on space 
or time. The boundary conditions at x" = —00 do not constrain further the 
coefficients 7^^", but they can be determined from the boundary conditions at 
a;° = +00. To achieve this end, we introduce the Fourier decomposition of the 
functions ci±f,xa(^)^ 

<fcAa(^) ^ E / {h^±icb(^";kXa) a%(2;) + /i^,(xO; kXa) a%^,ix)} . (38) 

It is then a simple exercise to rewrite the boundary conditions at .t" = +00 as 

E / b'-'^ f'^-LiPCb) 7I:'," h%^ipCb)] ^ 6{p q) /(+)(pC6) , 



e/K^." 



h^+LiPCb) - 7''," h^-LipCb)] = S{p - q) f^-\pCb) 



(39) 



where, to keep the expressions compact, we have omitted the argument x'^ = 
+00 in all the Fourier coefficients. This can be seen as a system of linear equa- 
tions for the coefficients 7±'^q . The solution of this system of linear equations is 
obtained in appendix A. 

Inserting the results in eq. (66) for 7^^" and Fourier decomposition in eq. (38) 
in eq. (37), one can easily determine the Fourier coefficients of b^. (x) at x'^ = 
+00 (eq. (35). Inserting these into eq. (30), we obtain^^ 



SHnJ^[z] 



EEE 



2 (2tt)^4:E„E„ 

LO;z(p),z(q) = l ^ ' P 1 \,a ^.b C,C 1^ 



5z{p)5z{q) 
X I {h^+LiP^h^-Li^Cc) + hi-LiP^b)h^+LiQCc)) f^+\pCb)f(+\qCc) 

+ {h%aipmi^LiQCc) + hi+^,Mb)h%M<^)) f^-Hpm^-HciCc) 
+ {hi-^xaipm'--LiQCc) + h'-_-Lipm%aiciCc)) f^+\pm^-\cicc) 

+ {h%aipm^-LiQCc) + hi^Lipmi-LiQC^)) f^-Hp(b)f(+\qCc)\ 

-(^^)^'5(p-q)^/(+)(pCc)/(-)(pCc) . (40) 

We have therefore obtained an expression for the connected piece of the two 
gluon spectrum entirely in terms of Fourier modes of the classical field (f^ ') 
and the small fluctuation field (^±feAa)- "^^"^ former can be determined by 



'Wc additionally use cqs. (65) to symmetrize the formula with respect to (p,q). 
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solving the Yang-Mills equations with retarded boundary conditions while the 
latter can be determined by solving the equations for small fluctuations about 
the classical field, also with retarded boundary conditions. 

The last term in eq. (40), proportional to 5{p — q) times the single particle 
spectrum, arises because the quantity dN2/d^pd^q is defined in such a way 
that its integral over p and q gives the average valuc"'^^ of N{N — 1). This 
term provides the —TV contribution to this quantity. Because the logs in the 
multiplicity N arise only at the order 0{g'^), this term cannot provide any 
leading log in the 2-gluon spectrum and can thus be dropped. 

3.4 Leading log resummation of the 2 gluon spectrum 

Combining the results in cq. (40) and eq.(28) in eq. (27), wc now have a formula 
for the 2-gluon spectrum, including both LO and NLO contributions. As men- 
tioned previously, it can in principle be evaluated, in full generality, by numerical 
solutions of small fluctuation partial differential equations with retarded bound- 
ary conditions. However, if one is interested primarily in the leading logarithmic 
piece of the NLO contributions, we can go signiflcantly further analytically. In- 
deed, as we will now show by using the information obtained thus far, we can 
compute the leading logarithmic contributions to the two gluon spectrum in 
perturbation theory. 

The first step in this derivation is to obtain an even more compact form 
for eq. (40) by using the linear operator T„ that we used previously in the 
expression for the 1-loop corrections to the single particle spectrum-see eq. (28). 
In Paper I, we demonstrated explicitly that this operator allows one to express 
the value of a retarded fluctuation at a point x in terms of the value of the 
classical field at the same point as 

a^'{x) = I d^u [a ■ T„] A^'ix) , (41) 

s 

where S is the initial surface on which we know the value of the fiuctuation. (The 
point X is located above this surface.) Performing the Fourier decomposition of 
both sides of this relation, we obtain simply the relation between the Fourier 
coefficients (at a;" = +oo) of the small fluctuation and the classical field to be 

h'-'\+oo;pXa) = j d^u [a • T„] /(')(+oo;pAa) . (42) 

s 

Applying eq. (42) to the various fluctuations that appear in eq. (40), and using 
the z{p) ~ 1 simplification of eq. (23), 



dN 
d^p 



SlnT[z] 



Sz{p) 



1 \-/(+)(pCc)/(-)(pCc), (43) 



(2^)32^;^ f- 



^^This can easily be checked on a Poisson distribution, for which the second derivative 
5\nJ^[z]/5z{p)Sz(q) is exactly zero. When we insert this in eq. (24) and integrate over p and 
q, wc obtain {N{N — 1)) = (A'')^ - as expected for a Poisson distribution. 
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it is a matter of simple algebra to check that 



5z{p)5z{q) 



-Sip-q) -jT- 



1.0;z(j,),z(g) = l 



dN 
(Pp 



C2 



dN 

connected d "p 



dN 
d^q 



(44) 



The subscript "connected" indicates that one of the T operators in the expres- 
sion £2 appearing in eq. (28) must act on the p-depcndent factor and the other 
on the q-dependent factor. (Terms where they both act on the same factor 
should be excluded.) 

We see now that eqs. (29) and (44) can be combined very easily, because the 
sum of "disconnected" and "connected" terms is equivalent to the unrestricted 
action of T„T.u on the product {dN/d^p){dN/d^q). We obtain thus 



d'^N. 



d^pd^q 










dN 


dN 


+ 


C1+C2 






LO 




d-^p 


d^q 



(45) 



where £1 and £2 were both introduced previously in eq. (28). 

We shall now discuss the logarithmic singularities in this expression. Firstly, 
6{p — q){dN/d^p)^^ does not contain large logarithms in x because these logs 
start appearing at NLO in the single gluon spectrum. Because we are restricting 
our discussion to leading logs, we can therefore discard this term henceforth. 
The logarithmic divergences in the second and third terms of the r.h.s. of 
eq. (45) can be extracted straightforwardly by using the main result of Paper I, 



^i+^^Lrog^H^^ 



Wi+ln 






Uo 



(46) 



Here ?ii,2 are the JIMWLK Hamiltonians of the nuclei moving in the +z and 
—z directions respectively [1,5-7], A represent the longitudinal momenta that 
separate the static color sources p\^2 in each of the nuclei respectively from 
the the gauge fields that produce gluons at the rapidity of interest, and M^ 
corresponds to the typical longitudinal momentum scales of the object (the two 
gluon spectrum in this case) to which the operator is applied. From eq. (46) we 
obtain 



d^iVo 



d^pd'^q 



LLoe 



1 + ln 



A+ 
M+ 



Wi+ln 



A- 



W2 



dN 
d^p 



dN 
d^p 



(47) 

All of our discussion thus far has been for a fixed distribution of sources pi,2 
in the two nuclei. The CGC effective theory [2-4,9-16], prescribes to average 
physical quantities over all the possible configurations pi2 of the fast color 
sources representing the projectiles, with gauge invariant weight functionals 
W^[pi.2] that describe the probability of each configuration. When we integrate 
eq. (47) over pi.2, we can exploit the hermiticity of the JIMWLK Hamiltonians 
Til. 2 in order to integrate by parts, so that the Hamiltonians are now acting on 
the distributions VK[/5i.2]. By reproducing the arguments developed in Paper 
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I for the single gluon spectrum, we obtain finally the factorization formula for 
inclusive two-gluon production, 



|^)^^_^/[d„][d«]»v.[..1v.-.,n|^ 



dN 



(48) 



at leading log accuracy. Here the distributions Ty[pi.2] obey the JIMWLK 
equation 

'-^^nwAp], (49) 

and are evolved thus from non-perturbative initial conditions at the beam ra- 
pidities to the rapidities Yi = ln(Y^/M+) and Y2 = \n{y/s/M~) respectively. 
In the regime where gluon radiation between the two tagged gluons is small, this 
formula resums all leading logarithms of l/xi.2 as well as all the rescattering 
corrections in (5^1,2)" to all orders. 

We now address the primary limitation of the present calculation. As the 
previous discussion hints, it is valid when the momenta p and q of the two ob- 
served gluons are close enough in rapidity so that they have similar longitudinal 
components. More precisely, we need to have 

a,ln(^j<l , a,ln^^j<l. (50) 

If this is the case, we can simply take M^ to be the common value-'^^ of p^,q^ . 
Physically, the condition of eq. (50) means that the probability of radiating a 
gluon between the two measured gluons is small. When the rapidity separation 
between the two gluons is large such that eq. (50) is violated, we need to resum 
gluon emissions between the tagged gluons; this would require a generalization 
of the present formalism, which is not discussed here. 

3.5 Factorization and the ridge in AA collisions 

A striking "ridge" structure has been revealed in studies of the near side spec- 
trum of correlated pairs of hadrons by the STAR collaboration [20-22]. The 
spectrum of correlated pairs on the near side of the detector (defined by an ac- 
companying unquenched jet spectrum) extends across the entire detector accep- 
tance in pseudo-rapidity of order A77 ~ 2 units but is strongly coUimated for az- 
imuthal angles Acj). Preliminary analyses of measurements by the PHENIX [23] 
and PHOBOS [24] collaborations appear to corroborate the STAR results. In 
the latter case, with a high momentum trigger, the ridge is observed to span 
the wider PHOBOS acceptance in pseudo-rapidity of A77 ~ 4 units. 

In Ref. [19], it was argued that the ridge is formed as a consequence of 
both long range rapidity correlations that are generic in hadronic and nuclear 
collisions at high energies, plus the radial flow of the hot partonic matter that 



^^It is of course not necessary that p+ and (jr+ be equal, just that they are close enough so 
that it docs not matter which value we chose between p+ and (jr+ . 
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is specific to high energy nuclear collisions. Let us first focus on the long range 
correlations that are essential to this picture-how are they generated? 

In the leading order formalism of the CGC, classical solutions of Yang- Mills 
equations are boost invariant [26-29]. Real time numerical simulations [SO- 
BS] also demonstrate that the Yang-Mills fields form flux tubes of a typical 
transverse size 1/Qs (where Qs is the saturation scale) with parallel chromo- 
electric and chromo-magnetic field strengths. (An important consequence is 
that these Glasma fields [25] have non-trivial topological charge [39].) Now, 
in section 3.1, we showed that the leading order 2-gluon spectrum, for a fixed 
configuration of sources, was given by eq. (26). Because each of the single 
particle distributions is boost invariant, the two particle spectrum is also, at 
this order, independent of the rapidity separation of the gluons. While the two 
gluons are uncorrelated for a fixed configuration of sources, correlations are built 
in through the averaging over the source distributions. In Ref. [19], the source 
distribution was assumed to be Gaussian as in the McLerran-Venugopalan (MV) 
model [2-4]. The ridge spectrum was shown to have the simple form 

/ dN\ 

-^^ = Cir, a) \'^W _ K, 

V7^ - ^^P'"^' I dN \ I dN \ ' a,(0,) ' ^^'^ 

\ dyp pdp d(t>p I \ dyq qdq d(f>q / 

where 

and Kiq is a number of order unity. For further details, we refer the reader to 
Ref. [19]. 

There are several conceptual issues in this context. Firstly, how does one 
justify this averaging procedure for the 2-gluon spectrum from first principles? 
Secondly, how docs one build in energy evolution of the sources? And finally, 
do NLO contributions spoil this picture? Our results in this paper solve most of 
these conceptual issues. Our result, in eq. (48), shows that the trivial LO result 
of eq. (26) can be promoted to a full Leading Log result simply by averaging 
it over the sources pi^2~with distributions of sources that evolve according to 
the JIMWLK equation. Most importantly, this shows that all higher order 
corrections, to leading logs in a:;i.2, do not spoil the form in eq. (51) of the 
Glasma flux tube picture and provides compelling evidence that it is a robust 
result beyond LO. As discussed previously, this picture will have to be modified 
when the rapidity separation between the gluons is greater than a^^((5s). 

These initial state considerations are not affected by the final state transverse 
flow of the Glasma flux tubes which is the other important feature determining 
the near side ridge seen in heavy ion collisions. Its been shown very recently 
that a proper treatment of flow and hadronization effects of the Glasma flux 
tubes provides excellent quantitative agreement with the RHIC data on the 
dependence of the ridge amplitude on ccntrality and as a function of energy, as 
well as the angular width of the ridge as a function of ccntrality [40] . Further 



16 



sophisticated treatments of both the initial state effects discussed here and the 
final state effects discussed in Ref. [40] , therefore open the door to quantitative 
3-D imaging of heavy ion collisions. A deeper relation between initial and final 
state effects, as outlined in Paper I, can be obtained by studying quantum 
fiuctuations at NLO, that are not accompanied by logs in a;i_2, but grow rapidly 
in time [41,42] in a manner analogous to plasma instabilities [43]. 

We should also mention that the initial state effects described here are also 
present in proton/deuteron-nucleus collisions [44-46], without the final state 
effects characteristic of the ridge in nucleus-nucleus collisions. These collisions 
are therefore useful in order to isolate the initial state effects and to corroborate 
the framework of multiparticle production in high energy QCD developed here. 

4 Multigluon inclusive spectrum 

In this section, we will show how the results of the previous section modify 
multi-gluon probability distributions, with the caveat, as previously, that these 
gluons are emitted in a narrow rapidity window. Wc will also derive a simple 
expression for the differential probability of producing n gluons. 

4.1 n-gluon spectrum at LO and NLO 

Our starting point in evaluating the inclusive n-gluon spectrum is eq. (8). Be- 
cause we have thus far obtained expressions up to NLO for the first and second 
derivatives of lnJF[z], it is convenient to rewrite this expression as^^ 

d"iV„ -A-(51nJ^[z] 






Because lnjr[z] = 0{g~'^) in our power counting, the LO term in the r.h.s. is 
of order g^^", the NLO term is of order g~'^'^^~^\ and NNLO and higher terms 
represented by the ellipses arc omitted at the level of the present discussion. 
The n-gluon spectra on the l.h.s. of eq. (53) are quantities that, for n > 1, are 
given by the first term on the r.h.s. By computing them to NLO we gain access 
to the first correction to the Poisson distribution, the deviation of the variance 
of the multiplicity distribution from the Poissonian result {N{N — 1)) = (A'')^ 
and the corresponding modifications for the higher moments of the distribution. 
We refer to the appendix C for a more detailed discussion of the interpretation 
of our result for the probability distribution of the gluon multiplicity. 



^This formula is obtained by replacing T[z\ by exp(lnjF[z]) in eq. (8). 
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At leading order, only the first term contributes, and wc obtain (for a fixed 
distribution of sources) 



d"A^n 



d3pi • • • d^p„ 
At next to leading order, we have 



n 



dN 



LO i = l 



d^P, 



(54) 



d''N„. 



d^Pi ■ ■ ■ d^p„ 



E 



dN 



(PPi 



n 

E 



dN 



d^p^ 
5^\nT[z\ 



Sz{p^)Sz{pA 



n 



dN 



d^Pk 



(55) 



All the objects that appear in this equation are known already from the discus- 
sion in Paper I and the previous section. In [1], we showed that 



dN 
d^p 



A+A 



dN 
d^p 



+ AiV„,o(p), 



(56) 



where Ci and C2 are defined in eq. (28). In the previous section, we showed 
that^"^ 



5'^\nJ^[z\ 



Co 



dN 

connected a J) 



dN 



(57) 



5z{p)Sz{q) ^^ 

where we remind the reader that the subscript "connected" attached to the 
operator £2 indicates that the two operators T it contains do not simultaneously 
act on the same object. 



4.2 Leading Log resummation 

If we combine the terms in eqs. (56) and (57), we get simply 



d"iV„ 



d3pi • • • d^p„ 



LLog 



Ci + £2] n 



dN 



(PPi 



(58) 



Using again eq. (46) and following the steps that lead from eq. (47) to 
eq. (48), we arrive at the all order leading log n-gluon spectrum 



d"A^„ 



d^p^ ■ ■ ■ d^p„ 



= J [Dpi] [DP2] Wy, [pi] Wy, [p2] 



dN 



(PPi 



dN 



d^p„ 



.(59) 



Once again, one needs all the rapidity differences between the n measured gluons 
to be much smaller than aj^, to ensure all leading logarithmic contributions 
are resummed by this formula. 



■''^We arc ignoring the term — 5(p — q) 



dp 



because it docs not contribute in the leading 



logarithmic approximation in x as discussed previously. 
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4.3 Generating functional in a small rapidity slice 

Eq. (59) provides a complete description of gluon production in the leading log x 
approximation when one considers a slice in rapidity of width Ay ^ aj^- One 
can summarize these results into a generating functional J^y,&y['^{p)] defined 
from the "master" JF[z(p)] as 



^y,^y[z{p)\ = Hz* {p)] with 




zip) if yp^[y-^y+^] 

1 otherwise 



(60) 

Setting the argument of the generating functional to unity outside of the phase 
space region of interest means that we define observables that are completely 
inclusive with respect to this unobserved part of the phase space. 

We see from eq. (9) that JF[z*(p)] can be obtained by multiplying eq. (59) 
by {z*{pi) — 1) • • • (z*(p„) — l)/n!, integrating over the n-gluon phase space and 
summing over n. Because z*{p) is unity outside of the strip of width AF in 
rapidity, the n-gluon spectrum outside of the regime of validity of eq. (59) is 
not needed. This procedure leads to a simple exponentiation of the leading log 
factorized formula for the generating functional T^ ^^ as 

(•^.,A.[^(P)]),,„^ = / [Dpi] [Dp2] W^y._-y [pi] W^y._+y [pa] 



X exp 



Y+i 



dPp (z{p) - 1) 



dN 
d^p 



(61) 



This leading log result for the generating functional, in turn, allows us to extract 
the corresponding formula for the differential probability of producing exactly 
n gluons in the rapidity slice of interest. This gives 



d"P„ 



d'^p^ ■ ■ ■ d^p„ 



= / [Dpi] [DP2] Wy,^..^^^y[pi] W^n_+r[p2] 



1 dN 



n! d^Pi 



dN 



d^Pn 



exp 



j3 dN 



Y-i 



(62) 



This simple result, valid, we emphasize, in the leading log approximation, sug- 
gests that the particle distribution in a small rapidity slice can be written as the 
average over pi_2 of a Poisson distribution with the leading log corrections com- 
pletely factorized into the JIMWLK evolution of the sources. Note that, despite 
appearances, eq. (62) is not a Poisson distribution after the integration over the 
sources, because particles produced uncorrelated in each configuration of pi and 
P2 are correlated in the averaged distribution because of the correlations among 
the color sources^*. 



^®For instance, two color sources may be correlated because they result from the splitting 
of a common "ancestor" in the course of JIMWLK evolution. 
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In general, even for a fixed distribution of sources, the probability distribu- 
tion is not Poissonian [17]. To some extent, the fact that we get a Poissonian 
functional form in the integrand of eq. (62) is a consequence of the way we 
have organized our calculation. In eq. (53) we are performing a weak coupling 

expansion of the moments {N{N — 1) {N — n + 1)), that includes the 

orders g~^" and the leading log part of the order g~^^"-~^\ Terms starting at 
the order g^2(n-2) ^^^ beyond the accuracy of our calculation, and therefore 
their value in our formulas are arbitrary. The arbitrariness of these subleading 
terms influences the precise form of the resulting generating functional. For 
example, if we had performed the weak coupling expansion of {N"') instead of 

{N{N — 1) {N — n + 1)), we would have obtained a different generating 

functional. Of course, the two generating functionals so obtained would lead to 
the same moments of the distribution to the order of our calculation. The non- 
trivial aspect of our result in eq. (62) is that all the deviations from a Poisson 
distribution that result from the large logarithms of x at NLO can be factorized 
into the JIMWLK evolution of the sources. Equation (62) shows how these cor- 
rections modify the n gluon production probabilities. The Poissonian nature of 
the multiplicity distribution and deviations from it are discussed in more detail 
in appendix C. 

5 Conclusion and outlook 

We demonstrated in this paper that our result of Paper I on initial state 
JIMWLK factorization for the single inclusive gluon spectrum in nucleus-nucleus 
collisions can be extended to inclusive multigluon spectra. Our result is valid 
provided all the gluons are produced in a rapidity window of width Ay < ctj^- 
Our final result for the generating functional for multigluon production, in the 
leading logarithmic approximation in xi^2j is very simple; the distribution of 
gluons produced in the stated rapidity window can be written as the average 
over the JIMWLK-evolved distributions of sources of a Poisson distribution. It 
is important to keep in mind that the result of this source average is not a 
Poisson distribution, due to the correlations between the evolved color sources. 

As we discussed in section 3.5, our results are of great interest in detailed 
imaging of the space-time evolution of nucleus-nucleus collisions. An important 
ingredient in future studies will be to extend the present result to the case of 
correlations between gluons produced at rapidity differences aj^ < AY. A full 
leading log computation of these initial long range rapidity correlations requires 
that one identifies and resums the additional large logarithmic corrections that 
may arise when the rapidities in the two-gluon spectrum are widely separated. 

An important caveat (also applicable to our previous study of the single 
gluon spectrum in nucleus-nucleus collisions) is that final state effects, related 
to the growth of unstable fluctuations, need to be resummed. While the de- 
tails are still unknown, the structure of the result is known. The result of the 
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resummation of unstable fluctuations, as shown in Paper I, can be expressed as 

(0>LLo8+u„.. = / [OAi] [DA^] H'V, [At] H',, [ij] 

x/[Z,aM]2|„M|0„,[JJ+„.A-+al^ (63) 

Here, we have traded the sources pi^2 in covariant gauge for the corresponding 
gauge fields Af2 — WJ-Pi.2- The functional Z[a(M)] is the spectrum of small 
fluctuations of the classical field on the forward light-cone. In Paper I, O cor- 
responded to the single inclusive spectrum but this formula also applies to the 
multigluon spectrum because the proof does not depend on the nature of the 
observable being measured. However, the complete functional form of Z[a{u)] 
is still unknown-for a first attempt, see Ref. [47]. 

These considerations are eased somewhat if we take the "dilute-dense" limit 
of proton/deuteron-nucleus collisions because we don't expect instabilities to 
play a major role in that case. Several studies have been performed in this 
limit [44-46,48-50]. A particular focus is on the applicability of the so called 
AGK cutting rules [51-54]. We plan to address these issues in a future work. 
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A Fourier coefficients of small fluctuation flelds 

We will outline here the solution to the system of equations 






,(+) f^/-h\_^Mahi+) 



":« h'^ilxaiPCb) 7^,," /»VfeAa(Kfo) = S{P Q) f'^> (pCb) , 



h^+LiPCb) - 7""" h^-LiPCb)] = S{P q) f^-HpCb) 



(64) 



that was obtained in eq. (39). We had previously derived analogous equations 
in the case of a simpler scalar theory in [17]. However, in [17], we did not 
manage to solve these equations and suggested that one may have to solve them 
numerically. It turns out that one can in fact obtain an analytical solution of 
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the eqs. (39), thanks to the relations 



X,a 



,(-) 



,(+) 






-h'-kXaiPib)h'+kXaiQCc) - (2^)^<54c '^bc 2i?p<5(p - q) , 



(+) 



,(+) 



(+) 



,(+) 



/^VfeAa(Pe^)/^-";Aa('ZCc) - /iL"fe\a(Pe^)/^VfcAa (^Cc) 



(-) 



,(-) 



(-) 



,(-) 



/^VfeAa(PCfc)/irfe^Aa('7Cc) - h'->,^{pCb)h^->,JqCc) 



(65) 



These relations are the mathematical consequence of the unitary temporal evo- 
lution of small fluctuations on top of the classical field A{x). In particular, 
an orthonormal basis of solutions of eq. (33) remains orthonormal at any later 
time. A proof of these formulas is presented in appendix B. Thanks to these 
relations, it is easy to invert the system of equations (39), and one gets 



kXa 



kXa 



7^4^ E [h+Li^Cb) f+\qCb) + h%JqCb) /(-)(gC&)' 



(66) 



B Unitary evolution of small fluctuations 

Consider the partial differential equation 



iy\ xdb 



(n.5/ - d,,d^,)s^ 



dU{Ae 



dA,a.{x)dA'if,ix) 



a"" (a;) = , 



(67) 



where we have written explicitly all the color indices. We assume that the 
background color field in which the wave propagates is real. For a generic 
solution a{x) of this equation, define the following vectors : 



\a) 



ai'^ix) 
d'"'lx) 



{a\ = (a*^"(x) a*^°(x)) , 



(68) 



where the dot means a derivative with respect to time. Then, it is trivial to 
check that the following "scalar product" , 



(ai|cr2|a2) =i.g^^(5ah / d^x dl^''{x)a2 {x) - a^^"(a;)a2 (x) 



(69) 
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where a^ is the second Pauli matrix, is independent of time when a^ and a^ 
are two solutions of eq. (67). 

Then, if the a±k\a{x) are the retarded solutions of eq. (67) whose initial 
conditions at x^ -^ — oo are e^(A;)T°e^*'^'^, one can check explicitly that 

(a+feAa|'7'2|a+fc/A'a') = {^TlflEk^W' ^aa' S{k - k') , 

{a-k\a\(T2\a-ki\'a') = -{2'n:f2Ek 5\\' 5aa' 5{k - k') , 

{o,+k\a\fJ'2\o,-kl\' a') = (oi-fcAa | Cr2 | Q+fc/A'a' ) = . (70) 

(Since all these scalar products are time independent, it is sufficient to check 
these relations by calculating the integral in the r.h.s. of eq. (69) for the corre- 
sponding initial conditions.) 

Consider now a generic solution a^(a;) of eq. (67). Since the solutions 
'^±feAa(^) span the entire space of solutions, we can write 

\a) = J2j [7^^" |a-fcAa) +7f " \a+kxa)] , (71) 

A,a ^ 

where the coefficients -y^^'^ do not depend on time. By using the orthogonality 
relations obeyed by the vectors \a±k\a), one obtains 

7^^" = -(a_feAa|T2|a) , 7f" = (a+feAa|'T2|a). (72) 

Inserting these relations back into eq. (71), one gets the following identity, 

I |a+fcAa)(a+fcAa| - |a-feAQ)(a-fcAa| = 0-2.9^" (5'"= , (73) 



E 



which is valid over the space of solutions of eq. (67). (The Lorentz indices /x, i' 
and color indices b,c do not appear explicitly in the l.h.s., but are part of the 
definition of the vectors la) and (al-see eq. (68).) This relation is valid at 
all times, and is the expression of the fact that the unitary evolution of small 
fluctuations preserves the completeness of the set of states \a±kxa) 

Let us now introduce states |a±fcAa)' ^^^^ ^^^ ^^^ analogue of the states 
Q±feAa) in the vacuum (i.e. when the background field is zero). Naturally, they 
are just plane waves aj^/.^^^ = e^{k)T°-e^^^'^ that we have introduced in order 
to perform the Fourier decomposition of classical fields and small fluctuations. 
The Fourier coefflcients h±jl^^{pC,c) of the fluctuations a'^^^ can be obtained 
as : 

(These relations are valid only in the regions where the interactions are switched 
off, i.e. when x° -^ ±00. In the rest of the discussion, we are only interested 
in these Fourier coefficients in the limit a;° -^ +00.) By multiplying eq. (73) by 
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(a°p^^|cr2 on the left and by (T2\a1,^^^ on the right and using (/ie/fcAa(PCc))* = 
''^-e'fcAa(PCc), we obtain the foUowing relation among these Fourier coefficients : 

^'" " -h^:klaiP^b)h^+kXai<l':c)] = <5«, e (27r)3 5^^,5bc 2ii;p<5(p - q) , 

(75) 
which is nothing but a compact way of writing the four eqs. (65). 

C Poisson distribution 

At first sight, eq. (62) appears to be the average over the distributions of sources 
of a Poisson distribution. This seems to contradict a result we stressed in 
[17], that the distribution of multiplicities calculated in a fixed configuration 
of sources pi^2 is not a Poisson distribution. For the sake of the discussion in 
this appendix, let us introduce the generating function F{z) for the multiplicity 
distribution in the region of rapidity \Y — Ay/2, Y + Ay/2] . In the language of 
the present paper, it is obtained by using in eq. (61) a constant function z{p) 
whose value is equal to the number z. 

Consider first this generating function for a given configuration pi2 of the 
external color sources. In [17], F[z) was parameterized as^^ 

oo 

lnF(z) = ^6^(z'--l) , (76) 

and we had obtained the formula for the probability P„ of producing n particles 
in the portion of phase-space under consideration to be 

n 
p=l ^ riH |-rp=n 

In Ref. [17], we also showed that h^ is the sum of all the cut connected vacuum 
graphs, where exactly r internal lines are cut. Because hj. is a sum of connected 
graphs, it has a perturbative expansion that starts at the order 1/ g^ , 

6r = — © 1 © 5^ ® • • • (78) 

In particular, all the bj. have a priori the same order of magnitude. However, it 
is easy to see that eq. (77) is a Poisson distribution only in the exceptional case 
where^" 

5i ^ , br = for r > 2 . (79) 



l^Compared to the notations used in [17], wc absorb the factors of l/g^ into the definition 
of the numbers b^ . 

■^"Froni eq. (77) and the definition of the fer, we have F{z) = J^^ z"Pn. Then, it is immediate 
to check that In F(z) should be a polynomial of degree one in the case of a Poisson distribution. 
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Since for a generic field theory, the br for r > 2 have no reason to vanish or to 
be smaher than bi , the distribution of the multipUcities in a fixed configuration 
of sources is in general not a Poisson distribution. Moreover, since 62.3,-- are 
of the same order in g^ as 61 , the deviations from a Poisson distribution is an 
effect of order unity, not a subleading correction. 

In order to make the connection with the present paper easier, it is preferable 
to parameterize F{z) as 

00 
lnF{z) = Y,Ckiz-l)'' . (80) 

k=l 

(This series starts at the index k ~ 1, because F{1) = 0.) The numbers Ck are 
related to the numbers bk by 

^^ = E(J(-i)'"'"^'»' ^^-Efj^- (81) 

k—r r—k 

where the (^) are the binomial coefficients. The derivatives of InF(z) evaluated 
at z = 1 are best expressed in terms of the coefficients Ck as 



d'' In F{z) 



dz^ 



= fe!cfc . (82) 



Let us now rephrase our results in this language. The inclusive n-particle spec- 
trum is the nth derivative of F{z) at z = 1. These derivatives read 

F«(l) = ci, 

F(2)(l)^c? + 2c2, 

F(3)(l) = cJ + 6ciC2 + 6c3,--- (83) 

All the coefficients Ck are sums of connected vacuum graphs, and therefore start 
at the order l/g^, up to logarithms. At Leading Order, we thus keep only 

LO 

At this order of truncation, one can obviously get a Poisson distribution, since 
this approximation is compatible with C2 = C3 = • • • = 0, i.e. 62 = ^3 = • • • = 
0. However, the coefficients 62, 3.-- could have any value of order g^^ without 
affecting our Leading Order truncation. The arbitrary choice one is allowed to 
make for these subleading terms in general alters the Poissonian nature of the 
distribution. 

The actual paradox arises only at the Next to Leading Order. There, one 
keeps the terms 

^^"^1)1 -h]"o+^NroHci]..o+^![ci]"-'N.o • (85) 

g-2(n-l) xlog 
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This does not correspond to a Poisson distribution anymore, since one needs a 
non-zero 62 in order to obtain these formulas. In fact, at this order of truncation, 
one has 62 = ci while the higher fe^'s are still zero. Even worse, our calculation 
of the second derivative of In F shows that ci is enhanced by a large logarithm, 
and is actually of order (/"^ ln(l/a;i_2) rather than the naive expectation g~'^ . 
Therefore, not only the distribution is not Poissonian, but the deviations from 
a Poisson distribution are logarithmically large. 

However, the main result of the present paper is that one can obtain the 
NLO corrections to the inclusive n-particle spectra by the action of a certain 
operator on the product of n 1-particle spectra at LO. In the present language, 
this reads 

^^"Hi) =[i + /:i + /:2]N^o. (86) 

LO+NLO 

Remember that so far all the discussion is for a fixed configuration of the sources 
pi, 2- Then, by averaging over these sources and by using the hermiticity of the 
operator £1 + £2, one can transfer the action of this operator from the quantity 
[ci]" to the distribution of sources. As wc have seen, this amounts to letting 
the distribution of sources evolve according to the JIMWLK equation. In other 
words, eq. (85) deviates strongly from a Poisson distribution, but does so in 
such a way that all correlations can be interpreted as coming from correlations 
among the sources that are generated by the JIMWLK evolution. 

Let us end this appendix with a word of caution in the interpretation of 
eq. (61). Strictly speaking, our Leading Log approximation gives us control 
only over the g^^ ln(l/a;i,2) part of the coefficient 62, but not over its g^^ part 
(without a log). The latter would only show up in a Next to Leading Log 
calculation. This means that in principle one could modify the argument of the 
exponential in the integrand of eq. (61) by a term of second degree in z(jp) — 1 
and with a coefficient of order g^^, without affecting any of our results for the 
inclusive gluon spectra at the order at which we calculate them. Obviously, 
such a modification of the integrand in eq. (61) would be a deviation from a 
Poisson distribution. Thus, the statement according to which the deviations 
from Poisson come from the JIMWLK evolution of the distributions of the 
sources pi.2 is true only for the largest of these deviations-i.e. those that are 
enhanced by large logarithms of the momentum fractions x\,2- Other deviations 
from Poisson exist, that are not enhanced by such logarithms-these are beyond 
the scope of the present calculation. 

References 

[1] p. Gelis, T. Lappi, R. Vcnugopalan, arXiv:0804.2630 [hcp-ph]. 

[2] L.D. McLerran, R. Vcnugopalan, Phys. Rev. D 49, 2233 (1994). 

[3] L.D. McLerran, R. Vcnugopalan, Phys. Rev. D 49, 3352 (1994). 

[4] L.D. McLerran, R. Vcnugopalan, Phys. Rev. D 50, 2225 (1994). 



26 



[5] L.D. McLerran, Lectures given at the 40'th Schladming Winter School: 
Dense Matter, March 3-10 2001, hcp-ph/0104285. 

[6] E. lancu, A. Leonidov, L.D. McLerran, Lectures given at Cargese Summer 
School on QCD Perspectives on Hot and Dense Matter, Cargese, France, 
6-18 Aug 2001, hep-ph/0202270. 

[7] E. lancu, R. Vcnugopalan, Quark Gluon Plasma 3, Eds. R.C. Hwa and 
X.N. Wang, World Scientific, hep-ph/0303204. 

[8] F. Gehs, T. Lappi, R. Vcnugopalan, Int. J. Mod. Phys. E 16, 2595 (2007). 

[9] J. Jalilian-Marian, A. Kovncr, L.D. McLerran, H. Weigert, Phys. Rev. D 
55, 5414 (1997). 

[10] J. Jalilian-Marian, A. Kovner, A. Leonidov, H. Weigert, Nucl. Phys. B 
504, 415 (1997). 

[11] J. Jalilian-Marian, A. Kovner, A. Leonidov, H. Weigert, Phys. Rev. D 59, 
014014 (1999). 

[12] J. Jalilian-Marian, A. Kovner, A. Leonidov, H. Weigert, Phys. Rev. D 59, 
034007 (1999). 

[13] J. Jalilian-Marian, A. Kovner, A. Leonidov, H. Weigert, Erratum. Phys. 
Rev. D 59, 099903 (1999). 

[14] E. lancu, A. Leonidov, L.D. McLerran, Nucl. Phys. A 692, 583 (2001). 

[15] E. lancu, A. Leonidov, L.D. McLerran, Phys. Lett. B 510, 133 (2001). 

[16] E. Ferreiro, E. lancu, A. Leonidov, L.D. McLerran, Nucl. Phys. A 703, 489 
(2002). 

[17] F. Gclis, R. Vcnugopalan, Nucl. Phys. A 776, 135 (2006). 

[18] F. Gelis, R. Vcnugopalan, Nucl. Phys. A 779, 177 (2006). 

[19] A. Dumitru, F. Gelis, L. McLerran, R. Vcnugopalan, arXiv:0804.3858 [hep- 
ph], to appear in Nucl. Phys. A. 

[20] J. Adams, et al., [STAR CoUaboration] Phys. Rev. Lett. 95, 152301 (2005). 

[21] F. Wang, [STAR Collaboration] talk at Quark Matter 2004, J. Phys. G 30, 
S1299 (2004). 

[22] J. Adams, et al., [STAR CoUaboration] Phys. Rev. C 73, 064907 (2006). 

[23] A. Adare, et al., [PHENIX Collaboration] arXiv:0801.4545 [nucl-cx]. 

[24] B. Wosiek, [PHOBOS CoUaboration], plenary talk at Quark Matter 2008, 
arXiv:0804.4352 [nucl-cx]. 



27 



[25] T. Lappi, L.D. McLcrran, Nucl. Phys. A 772, 200 (2006). 

[26] A. Kovner, L.D. McLerran, H. Weigert, Phys. Rev. D 52, 3809 (1995). 

[27] A. Kovner, L.D. McLerran, H. Weigert, Phys. Rev. D 52, 6231 (1995). 

[28] Yu.V. Kovchegov, D.H. Rischke, Phys. Rev. C 56, 1084 (1997). 

[29] M. Gyulassy, L.D. McLerran, Phys. Rev. C 56, 2219 (1997). 

[30] A. Krasnitz, R. Venugopalan, Phys. Rev. Lett. 84, 4309 (2000). 

[31] A. Krasnitz, R. Venugopalan, Phys. Rev. Lett. 86, 1717 (2001). 

[32] A. Krasnitz, R. Venugopalan, Nucl. Phys. B 557, 237 (1999). 

[33] A. Krasnitz, Y. Nara, R. Venugopalan, Nucl. Phys. A 727, 427 (2003). 

[34] A. Krasnitz, Y. Nara, R. Venugopalan, Phys. Rev. Lett. 87, 192302 (2001). 

[35] A. Krasnitz, Y. Nara, R. Venugopalan, Phys. Lett. B 554, 21 (2003). 

[36] T. Lappi, Phys. Rev. C 67, 054903 (2003). 

[37] T. Lappi, Phys. Rev. C 70, 054905 (2004). 

[38] T. Lappi, Phys. Lett. B 643, 11 (2006). 

[39] D. Kharzeev, A. Krasnitz, R. Venugopalan, Phys. Lett. B 545, 298 (2002). 

[40] S. Gavin, L. McLerran, G. Moschelh, arXiv:0806.4718 [nucl-th]. 

[41] P. Romatschke, R. Venugopalan, Phys. Rev. Lett. 96, 062302 (2006). 

[42] P. Romatschke, R. Venugopalan, Eur. Phys. J. A 29, 71 (2006). 

[43] S. Mrowczynski, hep-ph/0511052. 

[44] R. Baier, A. Kovner, M. Nardi, U.A. Wiedemann, Phys. Rev. D 72, 094013 
(2005). 

[45] J. Jahhan-Marian, Y. Kovchegov, Phys. Rev. D 70, 114017 (2004), 
Erratum-ibid. D 71, 079901 (2005). 

[46] K. Fukushima, Y. Hidaka, arXiv:0806.2143 [hep-ph]. 

[47] K. Fukushima, F. Gelis, L. McLerran, Nucl. Phys. A 786, 107 (2007). 

[48] M.A. Braun, Eur. Phys. J. C 42, 169 (2005). 

[49] C. Marquet, Nucl. Phys. A 796, 41 (2007). 

[50] A. Kovner, M. Lublinsky, JHEP 0611, 083 (2006). 



28 



[51] V.A. Abramovsky, V.N. Gribov, O.V. Kancheli, Sov. J. Nucl. Phys. 18, 
308 (1974). 

[52] M. Ciafaloni, G. Marchcsini, Nucl. Phys. B 109, 261 (1976). 

[53] J. Bartels, M. Salvadore, G.P. Vacca, Eur. Phys. J. C 42, 53 (2005). 

[54] J. Bartels, M. Salvadore, G.P. Vacca, JHEP 0806, 032 (2008). 



29 



